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Abstract 

A phenomenological kinetic energy theory of buoyant multiphase 
plumes is constructed, being general enough to incorporate the disso- 
lution of the dispersive phase. We consider an axisymmetric plume, 
and model the dissolution by means of the Ranz-Marshall equation in 
which there occurs a mass transfer coefficient dependent on the plume 
properties. Our kinetic energy approach is moreover generalized so as 
to take into account variable slip velocities. 

The theoretical model is compared with various experiments, and 
satisfactory agreement is found. One central ingredient in the model is 
the turbulent correlation parameter, called /, playing a role analogous 
to the conventional entrainment coefficient a in the more traditional 
plume theories. We use experimental data to suggest a relationship 
between /, the initial gas flux at the source, and the depth of the gas 
release. This relationship is used to make predictions for five distinct 
case studies. 

Comparison with various experimental data shows that the kinetic 
energy approach built upon use of the parameter / in practice has the 
order of predictive power as the conventional entrainment-coefficient 
models. Moreover, an advantage of the present model is that its pre- 
dictions are very quickly worked out numerically. 

Finally, we give a sensitivity analysis of the kinetic energy ap- 
proach. It turns out that the model is relatively stable with respect 
to most of the input parameters. It is shown that the dissolution is of 
little influence on the dynamics of the dispersed phase as long as the 
dissolution is moderate. 
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1 Introduction 

The motivation for studying bubble plumes is evident, from the fact that 
these plumes are encountered in a variety of engineering problems. To men- 
tion a few, plumes have been used to damp sea waves in harbours (pneumatic 
breakwaters), to prevent surface ice formation in harbours, to mix stratified 
fluid layers, to rc-acratc lakes, and finally to protect installations from shock 
waves produced by underwater explosions (cf., for instance, Bhaumik, 2005). 

With increasing subsea activities plumes have acquired increased impor- 
tance from a risk assessment point of view. It becomes imperative to obtain 
knowledge about the implications of a rupture, or even a breakage, of a sub- 
sea pipeline. Thus in case of an underwater blowout of inflammable gas, one 
wishes answers to the following questions: 

• What is the concentration of gas at the free surface? 

• What is the extension of the area covered by gas at the free surface? 

• What is the rising time of the gas? 

The aim of the present work is to provide a theoretical framework from 
which answers to these questions can be given, in a quick and reasonably 
accurate way. Our approach is based on the kinetic energy method, proposed 
by one of the present authors some years ago in the case of air-bubble plumes 
(Brevik 1977, 1996,1999). The essential generalization of the present paper 
is that we allow for dissolution of the gas. Therewith a realistic answer to 
the first question above becomes in principle obtainable. Especially when 
the plume is arising from large depths, the gas dissolution is expected to be 
an important factor. 

Figure 1 sketches a plume coming from a gas release at depth D. It is 
useful to divide the plume into three different zones: The zone of fiow estab- 
lishment (ZFE) is characterized by high gas concentration and fiow condi- 
tions changing rapidly with the height z above the real source (we thus do 
not introduce a virtual source being displaced some distance downward from 
the real source). Then, there is a surface zone, characterized by an outward 
radial flow of entrained fluid, yielding a rapid increase of the effective width 
of the plume. The third zone, the one of main interest here, is the estab- 
lished flow zone characterized by the self-similarity property, implying that 
the flow has established some sort of self-governing equilibrium (Socolofsky 
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et al., 2002). This means that the time averaged cross-sectional profiles of 
the plume parameters maintain a near-Gaussian form. The plume param- 
eters are fully determined by a centerline value and a width, both being a 
function of the height z. The Gaussian form of the plume parameters has 
been verified in several experiments; cf., for instance, Kubasch (2001). 




Figure 1: Sketch of rising plume in an unstratified environment originat- 
ing from a depth D below the (still water) surface. The plume consists of 
two boundaries, one containing the bubbles intermixed with entrained water 
(bubbly core, sketched with dashed lines) having a characteristic width Aa, 
and one denoting the outer boundary of upwards movement of water (water 
flow boundary, solid line) having a characteristic width a. The effect a ra- 
dial jet near the surface is sketched by an rapid increase in plume the plume 
width within the surface zone. Figure inspired by similar graphics found in 
Kubasch (2001) 
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As shown in Fig. [T] the radius of the plume is a, describing the water 
boundary, and the centerhne vertical velocity Uc- Furthermore, the dispersive 
phase is bounded by a characteristic radius Xa, where A < 1. The parameter 
A is introduced in order to account for the fact that the dispersed phase forms 
a "core", whereas the entrained fluid forms a wider conical structure. A is 
thus the relative width of the inner core compared to the outer water flow 
boundary. 

In the present paper we will be concerned with the zone of established 
flow only, as this zone is the dominant one in cases of moderate or large 
depths. We will ignore the phenomenon of detrainment (or peeling), which 
may occur in the case of stratification, or if the density of entrained fluid is 
significantly increased due to dissolution of the dispersive phase (i.e., the gas 
bubbles) . 

It is useful to distinguish between single- and multiphase plumes. The 
main difference being in the discrete nature of the buoyant dispersive phase. 
In a single phase plume the buoyancy is well mixed with the bulk fluid and 
the advection of buoyancy is controlled by the motion of fluid alone. For 
the multiphase plume the dispersed phase itself supplies the buoyancy to 
the plume and the distribution in the plume is controlled both by its own 
dynamics and by the motion of the bulk plume fluid. This distinction is es- 
pecially important in the case of stratification, and in the case of cross-flows. 
If cross-flows are present a separation between the dispersed and entrained 
phases can occur, yielding a behaviour qualitatively different from the char- 
acteristic plume behaviour. An extensive experimental study of multiphase 
plumes has been carried out by Socolofsky (2001). 

A great deal of effort has been laid down to describe the behaviour of 
multiphase plumes. The main focus of modelling has been on integral mod- 
els, considering the integral of relevant plume parameters over appropriate 
control volumes. The dynamics of the plume parameters is governed by cou- 
pled differential equations derived from first principles and/or closure models 
for turbulence, mass transfer, heat transfer, and entrainment. 

According to the entrainment hypothesis, as introduced by Morton et 
al. (1956), the rate of entrainment at the edge of the plume is proportional 
to some characteristic velocity at that height. The characteristic velocity is 
taken to be the centerline velocity Uc at that height, and the proportionality 
constant is the entrainment coefficient a. Mathematically, 



dz 



2TT(jaUc, (1) 
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where Q^, is the volumetric flux of cntrainmcnt water. 

Two different strategies are followed when dealing with multiphase plumes. 
They depend in turn on the way in which the multiphase nature of the prob- 
lem is treated. A single plume model (or mixed-fluid model) treats the plume 
as a mixture of dispersed phase and entrained fluid. If the ambient fluid can 
be considered as a stagnant and unstratisfled fluid, one calls this a simple 
plume. 

Simple plume models have been considered in various contexts. Morton 
et al. (1956) used a simple plume model in their paper introducing the 
entrainment hypothesis; Milgram (1983) compared his experimental results 
with the values predicted by a simple plume model, whereas Wiiest et al. 
(1982) used both a simple plume model with the entrainment hypothesis to 
describe the effect of dissolving plumes. The work of Brevik (1977) is also 
based on a simple plume formulation, though without making use of the 
entrainment hypothesis. 

A different strategy used to deal with the multiphase nature of the plume 
is to use the concept of two-fluid models, namely to treat each phase sep- 
arately and introduce coupling terms in the differential equations for each 
phase. This superimposition of two plume-like structure leads to a formalism 
called a double-plume m.o<\e\ (or two-fluid model). In the past, two- fluid mod- 
els have been implemented only together with the entrainment hypothesis. 
McDougall (1978) modelled the system as an inner circular plume containing 
the gas bubbles and some entrainment water, and an annular plume contain- 
ing only water. The formalism of McDougall was generalized by Asaeda 
and Imberger (1993), and by Crounse et al. (2007), simplifying the imple- 
mentation of density effects of dissolving bubbles. An extensive comparison 
between mixed and two-fluid models based on the entrainment hypothesis 
is given by Bhaumik (2005). The advantages of the integral models are 
that the governing equations allow insight into the flow dynamics; they are 
computationally efficient and produce reasonable results in many cases. A 
drawback is that the integral models lose their vahdity as the plume becomes 
less self-similar (Socolofsky et al, 2002). 

In the next section we summarize the key properties of the kinetic energy 
model. In section 3 we solve the nondimensional governing equations and de- 
termine the values of the parameters, especially the correlation parameter 7, 
so as to get a reasonably good agreement with existing experiments. Section 
4 is devoted to large-scale case studies, of importance for the oil industry. In 
section 5 we analyze the sensitivity of the introduced parameters. 
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2 Model formulation 



The model forming the basis of our present investigation is the kinetic energy 
model presented by Brevik and Killie (1996) for axisymmetric air-bubble 
plumes, now generahzed so as to account for dissolution of the dissolved 
phase. Here the semi-empirical Ranz-Marshall equation (1952) is important 
for predicting the gas dissolution. This is moreover similar to the approach of 
Wiiest et al. (1992); they describe dissolving plumes using the entrainment 
hypothesis. Our model is based on a mixed-fluid formalism, and intends to 
describe a simple bubble plume in the established zone. We shall predict the 
width, dissolution and the rise time of a plume originating from the source 
lying at depth D (^ = 0). 



2.1 Mass balance 

Consider first the gas dissolution. The rate of mass transferred from the 
dispersed phase to the ambient fluid follows from the Ranz-Marshall equation 
(1952): 

rim 1 

-47rrlK{cs-Ci), (2) 



dt 

where is the mass of one spherical bubble with radius Vf,, Cg is the solubility 
of the dispersed phase, q is the local concentration of the dissolved species 
in the ambient fluid, and K is the mass transfer coefficient. 

Whereas the original Ranz-Marshall equation was derived to describe the 
evaporation of pure liquid drops, the equation has later been successfully 
applied to describe dissolution of gas bubbles (Wiiest 1992, Crounse 2007, 
Johansen 2000), and has become a standard in this area of research. The 
coefficient K depends on bubble size and flow parameters, as discussed by 
Wiiest (1992), Zheng and Yapa (2002), and Crounse (2007). 

Zheng and Yapa (2002) combine equations from earlier works to derive a 
general formula for K for bubbles in water. Assuming spherical bubbles of 
radius r^, 

K = 0.0113./^^^ (3) 

where ro = 0.45 cm, Us is the bubble slip velocity, and V is the molecular 
diffusivity of gas in the liquid. 
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As the mass of a bubble is = (47r/3)r^pg with pg the gas density, we 
may insert this into ([2]) and moreover make use of the equation 



P,ir.) = (4) 



atm 



which follows from the assumption of isothermal expansion. Here Pgs is the 
gas density at the surface, Patm is the pressure at the free surface expressed 
as a head of water column, and D = Patm + D- This leads to the following 
equation for the height dependence of the bubble radius: 

dn ^ KPaimjCs - Cj) ^ Tb 

Uzb being the vertical velocity of a bubble. 

Consider next the mixed-fluid continuity equation by summing the stan- 
dard continuity equations for each phase over the two phases. The density p 
of the mixture is 

p = aipi + agPg, (6) 

where ai, ag are the phase fractions of the liquid and gas phases. Assuming 
that ag <^ ai (or ai ^ 1), the last term in is negligible. Furthermore, we 
assume 

p^ pi^ p^, (7) 

where p^ is the ambient fluid density, taken to be constant. Consequently the 
(stationary) continuity equation for an axisymmetric plume can be written 

as 

-^{rUr) + TrUz = 8 
r or oz 

2.2 Momentum equation 

Performing the same kind of averaging procedures as in earlier papers (Brevik 
and Killie 1996, Brevik and Kluge 1999), now including gas dissolution, and 
using the approximation ([7]) meaning that all momentum contributions from 
the dissolved phase are neglected apart from the buoyant term, we can write 
down the z component of the momentum equation for the fluid mixture. We 
assume that the plume is fully turbulent, so that the velocity field can be 
Reynolds decomposed into a mean value and a fluctuating term. Viscous 
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stresses are small and negligible compared to turbulent stresses. For the 
averaged Navier-Stokes equation we obtain, in standard notation, 



d d d d d 

Pw[ur—u^ + + Uz-^^u^ + = --Q^P - Pw9 + PwOigg (9) 

where g is the gravitational acceleration. The mean pressure p can be written 
as a sum of a dynamic and a hydrostatic part, p = Pd + Pwg{D — z) . Thus the 
right hand side of ([9]) simplifies to PwCtgg. Introducing the Reynolds stresses. 



Trz = -pwKK^ 'Tzz = -pwU'^, (10) 

and making use of the continuity equation ([H]), we can write the correlations 
in (ED as 



where it is assumed that the lateral variations of the Reynolds stresses are 
dominating. 

The momentum equation thus reduces to the form 

do 1 d , _ . ^ d . , 

—u^ + - — {rUrU;,) = ttgg H ^(rT:,J. (12) 

oz r or pwT or 

This equation is integrated over a cylindric control volume having the bound- 
aries z = [0,z'] [z' arbitrary), and r = [0, oo]. The boundary conditions are 
summarized in Tabled! 

Let Pn be the number density of bubbles. The void fraction ag is related 
to Pn via 

47r o , . 

With the given boundary conditions, the integration of (JT^D leads to 

roo rz' roo SvT^ C^' 

271 / u^rdr = 27rg / / agvdzdr = g / puvdzdr. (14) 

Jo Jo Jo 3 JO Jo 

In order to advance further, assumptions must be made about the form of 
Uz and Pn- We shall assume Gaussian distributions: 

Uz{z,r) = Ua{z)e 5^ (15) 
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Table 1: Boundary conditions for integration of momentum and kinetic en- 
ergy equations 



Limit 


Condition 


r = 


Ur = 




rTrz = 




Uz : finite 


r ^ oo 


Uz = 




rTrz 




rUr : finite 



2 



pn{z,r) = pnciz)e (16) 

where Uc{z) and Pnc{z) are (unspecified) centerlines values. The plume grows 
wider as the free surface is approached, mainly due to turbulent diffusion, 
which corresponds to cr = a{z). 

Insertion of (fT5l) and (fT6l) into f|T^ leads to the relation 



ula^ = —gX^ / r^^pn^dz. (17) 

6 Jo 



Upon differentiation. 



-^ula^ = ^gX'^rlpncCT^ = 2gX''VbPncCr'' (18) 

where Vb = Vb{z) is the volume of a single bubble. 

Assuming that the number of bubbles is constant, in other words that 
there is an equilibrium between bubble coalescence and break up, the center- 
line number density p„c can be expressed in terms of the number flux density 
N of bubbles. Evaluating the basic expression 

N = 2tt I pnUbrdr, (19) 

we obtain 







N{l + X^) 
" 27rAV2K + M,(l + A2)] 

Here the factor Uc + Us{l + X'^) may be interpreted as a typical vertical bubble 
velocity, i.e. Uzb = Uc + Us{l + A^). 
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Inserting (120]) into (ITS]) we obtain the final form 
d gV,N{l + X') 

It is desirable here to check with the known result in the case of no gas 
dissolution. Putting iiT = in ([5]) and integrating with respect to z, we get 

^=(^^)^ (22) 

Cubing this equation we get 

H = H.^. (23) 

D — z 

which after insertion into (12TI) yields 

d 2 2 gV,sNPatra{l + \^) qQ"^ Patm{l + >?) 
—UM = z = z (24) 

dz Tr{D-z)[u, + Us{l + X^)] 7r{D - z)[uc + Us{l + X^)] 

where = VbgN agrees with Eq. (22) in Brevik and Killie (1996). 

2.3 Kinetic energy equation 

Equations and (1^ contain three unknown quantities, Uc, T}, and a (as- 
suming the other quantities can be determined from second principles or 
experiments). In order to close the system of equations, a third equation is 
needed. In our approach, this equation will be the conservation equation for 
kinetic energy. 

The starting point is the momentum equation in the form (fT2l) . Taking 
again into account the continuity equation ([8]), we obtain after some manip- 
ulations 

d A 1 d ,1 _ Uz d , , _ 

This equation corresponds to (fT2l) in the momentum case. We integrate (125!) 
over the same control volume as before, and apply the boundary conditions 
given in Table 1. Then differentiating the equation with respect to z we get 
the integral-differential equation 

roo \ roo Q roc 

— / -ulrdr = / ru'u'—Uzdr + g / a„Uzrdr. (26) 
dz Jo 2 Jo or Jo 
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In addition to the models for Uz and a^, a model for the turbulent correlation 
u'j.u'z is needed. As in Brevik and Killie (1996) we assume self-preservation: 



fiv), (27) 



where / is an unspecified function of the nondimensional parameter rj = r/cr. 
making use of this, together with the Gaussian forms (fTSjl and (fT6|l . we obtain 
finally 

d 3 2 3 gVbU^N 3 
dz TrUc + Us[l + X) 

where 

J = 6 / /r/^e-a" rfr] (29) 





is an unspecified constant to be evaluated from experiments. 

It is easily shown that (128|) is consistent with the results of Brevik and 
Killie (1996) in the limit K ^ 0. 

Equations ([5]), (|2T1) and (1281) form a closed set for the unknowns r^, and 
a. For given initial conditions, if the parameters Ug, A and / are determined 
from experiments, the equations determine the unknown quantities at an 
arbitrary height within the zone of established flow. 



3 Method of solution 

3.1 Nondimensional formulation 

For computational purposes it is convenient to express the equation in nondi- 
mensional form. We define the parameters 

* ^ if K 

z = = — , K = — , 

D D Wq 

u: = ^, < = — , y = (30) 

^0 ' Wo pggD 

where Wq is a reference velocity defined in section 3.3. Equation ^ can now 
be written 

*1 1^ rl .3^^ 

dz* [ul + {I + \^)ul){l - z*) 3(1 -;z*) ^ ' 
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Further, and ([2HD become 



d 



(rtni + A^l 



dz*^ " ' ~ u* + u*Al + A2) 



dz* 



u* + u*{l + A2) 



with 



2 \ D^gN 



In ( [33|) . we have introduced 



cr 



Further introducing 



Scr^u^Wo 



AD^gN 

we may after some manipulations replace fl?Il) and 



with 



*\2 



+ + A2)) 



c/ 



(32) 
(33) 

(34) 



(35) 



(36) 



(37) 



(38) 



M* +M*(1 + A2)' 

The three nondimensional equations (13T!) , (1371) and (l38l) are solved iteratively 
using the ODE45 solver in Matlab. This solver is a Runge-Kutta method for 
numerical integration which uses variable time steps for efficient computation. 
We solve the problem as an initial value problem. 



3.2 Initial conditions 

The structure of the governing equations makes them sensitive to the initial 
values of the variables. Care should therefore be taken to determine these 
values accurately. In particular, as (!37|) is singular for n* = 0, finite initial 
conditions for Uc and a need to be identified. 
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Bhaumik (2005) compares three different concepts for determining the 
initial conditions: 

• Power series (McDougall, 1978) 

• Virtual point source (Ditmars and Cederwall, 1974) 

• Densimetric Froude number (Wiiest et al., 1992). 

The conclusion of Bhaumik is that the method of Wiiest et al. is superior 
to the others as it preserves the multiphase nature of the plume through the 
incorporation of the phase fraction Ug and relative width A and is indepen- 
dent of any arbitrary parameters. This method is accordingly used in the 
following. 

As for the critical conditions on r^, the main parameter controlling the 
bubble size is the pore size of the source (or diffuser) ( Wiiest 1992, Bhaumik 
2005). These authors point out that the the initial bubble size depends on 
the gas flow rate, larger flow rates yielding larger bubbles and vice versa. The 
exact correlation between flow rates and bubble sizes is however not known. 
For simplicity, we therefore assume in the following that the conditions are 
ideal in the sense that the initial bubble sizes are determined solely from the 
orifice diameter if not provided directly from observations. 

Next, consider the initial condition on a. We will follow Wiiest (1992) in 
taking the initial plume radius to be given by the apparent size of the bubble 
source area, thus considering the entire diffuser uniform source. The 

diffuser area equals the initial area covered by the core. 



where 6° is the initial width of the top-hat distribution used by Wiiest, and 
A is the relative with introduced earlier. In our notation, the top-hat width h 
is equal to the standard deviation a (this comes from the condition that the 
number and momentum fluxes of the plume should be independent of which 
distribution is used). Thus IP = 2o"° initially, and fl5^ yields 



where is the apparent radius of the diffusor system. 

Consider finally the critical centerline velocity Following again Wiiest 
(1992), we introduce first a densimetric Froude number Fr which in our 
notation reads 



(39) 



Fr 




(41) 
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As the density p can be expressed in terms of phase fractions, 



P = (^gPg + (1 - (^9)9-' 



(42) 



and as the void fraction ag can be related to the volumetric flux Qg of gas, 



mation, the initial condition on vP^ is determined by Fr°. 

3.3 Parameter determination 
Relative distribution width A 

The parameter A determines the relative width of the velocity distribution 
and the number density distribution. The number density profile is neces- 
sarily bounded by the velocity distribution profile (as this defines the outer 
boundary of the plume). The parameter A should thus take a value on the 
interval [0, 1]. 

The literature presents several values for A, ranging from A = 0.2 used 
by Ditmars and Cederwall (1974) to A = 1 used by Crounse (2007). Wiiest 
(1992) used the value A = 0.8. The value of 0.8 is also adopted in the 
comparative study of existing models conducted by Bhaumik (2005). The 
parameter A is taken to be 0.8 in the present work as the work by Wiiest 
(1992) and Bhaumik (2005) showed that it yields good results for dissolving 
plumes. 

Bubble slip velocity Ug 

For a particle moving with a steady terminal velocity C/ in a gravitational 

field, the drag force balances the difference between weight and buoyancy 
(Clift (2005)). It is reasonable to expect that the actual slip velocity for 
the bubbles is different from the terminal velocity because of screening and 



Qg = 47ra^ag[uc + Us{l + A^)], 



(43) 



we obtain, using ag <ti ai, that 
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interactions between bubbles as shown in simulations carried out by Esmaeeli 
and Tryggvason (1999). However, in the present work it is assumed for 
simplicity that the bubble slip velocity equals the terminal velocity. 

The curve fitting obtained by Wiiest (1992) for terminal velocities is used 
in the present work and is given by 

usin) = wi(^)^-35' if 7 < 7-0 • 10-^ (45) 

Usin) = Wo if 7.0 • 10-^ <j<5.1- 10-^ (46) 

Us{n) = W2Cjf' if ^> 5.1 -10-3 (47) 

where wi = 4474 m/s, wq, = 0.23 m/s, W2 = 4.202 m/s and f = 1 m. wo 
is the reference velocity introduced in the dimensionless formulation. 



Dissolution parameters 

The Ranz-Marshall equation introduces three new parameters K, Cg and q. 
As shown, the mass transfer coefficient K can be determined from bubble 
properties, given that the molecular diffusitivity V is known. The molecular 
diffusitivity for gases in water is of order 10~^ cm^/s (Perry and Green 1997) 
and is general dependent on viscosity (//) and temperature (T). As the details 
of these parameters would complicate the model further, the generic value of 
T> = 10~^ cm^/s is chosen. 

The solubility Cg is in general a complicated function of the thermody- 
namical state of the system. For moderate pressures (up to 50atm), the 
solubility can be expressed by Henry's law: 

Cs = Hp (48) 

where H is Henry's constant and p is the partial pressure of the phase in 
question. H is material and temperature dependent. Suitable values for H 
are found in comprehensive handbooks such as Lide and Predrikse (1995). 

The rate of mass-transfer from a bubble is dependent upon the concen- 
tration of dissolved species in the surrounding water q. The process of dis- 
solution tends to increase this concentration. Assuming that the increase in 
concentration is of the order of the initial concentration of the ambient fluid, 
and that this concentration is small compared to the saturation concentra- 
tion, the in-situ concentration of dissolved gas q is negligible compared to 
the solubility c^. 
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The partial pressure of the gaseous phase is for simphcity approximated 
with the hydrostatic pressure in the surroundings, thus neglecting effects of 
surface tension in the bubble. 

3.3.1 Turbulent correlation parameter I 

The turbulent correlation parameter / introduced is based on the mean cor- 
relations between turbulent velocity components, i.e. a Reynolds stress and 
is in the present work determined from experimental data. The experiments 
of Milgram (1983) are chosen for calibration. The experiments were carried 
out with air released in an unstratified environment with negligible effects of 
crossfiow and are thus similar to the assumptions used in the present model. 
Four different initial flow rates Qo were used, varying from 0.0312 kg/s to 
0.767 kg/s. 

Using the given input data and initial conditions, the model equations 
are solved with J as a free parameter. The parameter I is varied until good 
accordance with the experimental data for the centerline velocity Uc and 
plume width h is obtained (6 = \/2a). Sample plots are given in figures [2] 
and [3] and overall results are presented in tables [2] and [31 

Table 2: Key results from calibration from Milgram (1983). Four release rates 
Qo are investigated from a release depth D of 50 m. The table gives relevant 
quantities for simulations done with optimal values for I (lopu) chosen. The 
mean centerline velocity Uc is calculated from the depth of release {D=50 m) 
divided by the total rise time t^se- 



^surf '^surf ^rise ^opti 

(kg/s) (kg/s) m (s) (m/s) ( ) 

0.031 0.027 3.368 67.460 0.741 0.075 

0.153 0.137 4.598 47.894 1.044 0.102 

0.368 0.331 6.563 45.500 1.099 0.147 

0.767 0.703 6.729 35.973 1.390 0.151 



As figures [2] and [3] show, the model presented yields satisfactory results 
when compared to the experiments conducted by Milgram (1983), if the free 
parameter / is chosen correctly. The plots do however show some deviations 
close to the surface. The deviations are due to plume interaction with the 
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Plume width b (m) 

Figure 2: Determination of optimal / from b, Qo = 0.7670 kg/s. Solid line 
shows optimal fit for the given mass rate (based on optimum for b and Uc), 
while dashed lines shows results when I is chosen ±5, ±10 and ±20% from 
optimum. 
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♦ Milgram (19B3) 
1 = 0.1208 

1 = 0.1359 

1 = 0.1435 

1 = 0.1510 

1 = 0.1586 

1 = 0.1661 

1 = 0.1812 



°1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 
Centerline velocity (m/s) 

Figure 3: Determination of optimal / from Uc, Qo = 0.7670 kg/s. Solid line 
shows optimal fit for the given mass rate (based on optimum for b and Uc), 
while dashed lines shows results when I is chosen ±5, ±10 and ±20% from 
optimum. 
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Table 3: Key results from calibration from Fannel0p and Sj0en (1980) (Data 
reproduced by Milgram (1983)). Four release rates Qq are investigated from 
a release depth D of 10 m. The table gives relevant quantities for simulations 
done with optimal values for / (lopu) chosen. The mean centerline velocity 
Uc is calculated from the depth of release {D=10 m) divided by the total rise 
time trise- 



Qo 


Qsurf 


bsurf 


trise 


Uc 


I opti 


(kg/s) 


(kg/s) 


m 


(s) 


(m/s) 





0.007 


0.006 


1.002 


12.770 


0.783 


0.100 


0.013 


0.013 


1.193 


11.291 


0.886 


0.120 


0.019 


0.019 


1.241 


10.085 


0.991 


0.125 


0.029 


0.028 


1.385 


9.498 


1.053 


0.140 



surface, yielding a radial jet of entrained water, described for instance by 
Brevik and Kristiansen (2002). Effects in the surface zone are however not 
considered in the present work. 

Brevik and Killie (1996) point out that the turbulent correlation param- 
eter / in general could be a function of the initial mass flux Qq and depth 
of release D. Tables [2] and [3] show that the value of / increases with increas- 
ing Qo- However, even though the values of Qq in table [3] are significantly 
lower than in table O the optimal values of / are similar. This suggests that 
decreasing values of D yield increasing values of /. 

It is not possible to combine D and Qq into a dimensionless group. This 
suggests the presence of some other parameter also influencing the turbulent 
correlations. A natural choice is the viscosity coefficient /i. Even though 
viscosity is neglected in the governing equations, it is still important for 
damping out turbulence at small scales. The following dimensionless group 
can be constructed based on the three mentioned parameters, Qq and fi: 

Re.^9±. (49) 

The dimensionless group Rce behaves as a Reynolds number Re\ 

Qq puV 

Rce = — = (50) 

Dfi fj, 
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where V is some length scale given by the plume cross section. In order 
to make predictions for an arbitrary value of Rce-, the values presented in 
tables [2] and [3] are fitted to a function of the form / = alni^e^; + 6, using the 
least square scheme Isqnonlin in Matlab. The function is chosen because of its 
relatively simple form and its ability to qualitatively reproduce the behaviour 
observed. The purpose of this approach is to provide a quantitative formula 
for predicting / and to aid in achieving a qualitative physical understanding 
of how relevant parameters influence /. 

The results of the least square fit is presented in figure HI 

0.1S 
0.16 
0.14 
0.12 
0.1 
O.OS 
0.06 
0.04 
0.02 


2 4 6 S 10 12 14 16 IS 20 
ReE=Q„/(Dn) 

Figure 4: Least square fit for relation between / and Rce- Solid circles show 
values for / obtained from the experimental data of Milgram (1983) and 
stars show values for / obtained from experimental data of Fannel0p and 
Sj0en (1980). Solid line shows the least square fit to a function of the form 
J = a In Re e + b. 

Optimal fit is obtained by choosing coefficients a and b corresponding to: 
/ = 0.0219 In i?ei5 + 0.1010. (51) 
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4 Case studies 



In the following, a release of natural gas consisting of a mixture of 85% 
methane and 15% ethane (mass basis) will be simulated from release depths 
of 100 and 300 m. For both release depths mass rates of 10, 25 and 50 kg/s 
will be considered. The gas mixture is assumed to behave as an ideal gas 
and is assumed to expand isothermally. The sizes of bubbles will be varied 
between 1 and 10 mm and the ambient temperature and the temperature of 
the gas will be held constant at 278 K. The solubility of gas is calculated 
from Henry's law (equation HHl) and the flow is assumed to be critical at the 
source and the size of the source is thus determined from the given mass rate. 

The turbulent correlation parameter / is determined from the initial rate 
of release Qo and the depth of release D by equation [SH As mentioned, 
equation [5T] is meant as a first estimate, and experimental data exist only 
up to values of -Re^; ~ 20. However, the test cases that are interesting for 
the industry require extrapolation well beyond this value. To deal with this, 
a choice is made to not allow for values of Rce predicting values of I larger 
than 35% of the observed values for /. The case of a 50 kg/s release from 
100 m falls outside this limit, and is thus omitted from the analysis. Results 
from the case studies are presented in table HI 
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Table 4: Key results from case studies. Mass flux at surface Qsurf, plume 
width at surface bsurf and plume rise time t^ise are computed from theory 
and mean centerline velocity Uc is calculated from the depth of release D, 
divided by the total rise time t^ise- 



D 


Qo 


' b 


Qsurf 


bsurf 






(m) 


(kg/s) 


(mm) 


(kg/s) 


(m) 




(m/s) 


100 


10 


1.0 


4.61 


18.22 


51.18 


1.95 


100 


10 


3.0 


8.01 


17.48 


47.38 


2.11 


100 


10 


5.0 


8.71 


17.39 


46.83 


2.14 


100 


10 


7.0 


9.00 


17.36 


46.76 


2.14 


100 


10 


10.0 


9.23 


17.33 


46.76 


2.14 


100 


25 


1.0 


13.96 


19.75 


39.07 


2.56 


100 


25 


3.0 


20.99 


19.16 


36.85 


2.71 


100 


25 


5.0 


22.40 


19.09 


36.52 


2.74 


100 


25 


7.0 


22.95 


19.06 


36.48 


2.74 


100 


25 


10.0 


23.42 


19.04 


36.48 


2.74 


300 


10 


1.0 




0.00 






300 


10 


3.0 


2.01 


48.91 


315.02 


0.95 


300 


10 


5.0 


4.51 


45.82 


285.10 


1.05 


300 


10 


7.0 


5.64 


45.01 


278.11 


1.08 


300 


10 


10.0 


6.66 


44.44 


273.57 


1.10 


300 


25 


1.0 


0.00 


72.78 


404.86 


0.74 


300 


25 


3.0 


7.91 


52.38 


234.92 


1.28 


300 


25 


5.0 


13.56 


50.17 


218.66 


1.37 


300 


25 


7.0 


15.98 


49.52 


214.57 


1.40 


300 


25 


10.0 


18.11 


49.05 


211.90 


1.42 


300 


50 


1.0 


0.01 


75.62 


297.89 


1.01 


300 


50 


3.0 


20.14 


55.24 


189.20 


1.59 


300 


50 


5.0 


30.26 


53.50 


178.84 


1.68 


300 


50 


7.0 


34.50 


52.96 


176.19 


1.70 


300 


50 


10.0 


38.18 


52.56 


174.41 


1.72 
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5 Discussion and sensitivity analysis 



In order to investigate the robustness of the model, a sensitivity analysis of 
critical parameters is conducted. Sensitivity analyses for A, i?D, r°, Fr°, 
T> and / are performed keeping other variables constant. In the sensitivity 
analysis the effects of parameter variation on the plume width 6, mass-fiux 
at the surface Qsurf and the mean centerline velocity Uc are investigated. 

The prediction of the mass-fiux of the dispersed phase is relatively ro- 
bust. Some variation (~ 1%) is found when varying the bubble radius rj,, 
(smaller bubbles yielding more dissolution) and when varying the molecu- 
lar diffusitivity T> (larger diffusitivitics yielding more dissolution). These 
findings are in accordance to those found by Wiiest (1992). The width of 
the plume b and mean centerline velocity Uc show little sensitivity to these 
parameters, suggesting that the dissolution itself has little influence on the 
other plume properties, as long as the dissolution is moderate. 

The analysis shows that the parameters of interest are insensitive to vari- 
ations of the diffuser radius and the dcnsimctric Froudc number Fr'^. 
The Froude number Fr^ does however have some influence on the centerline 
velocity Uc within the zone of flow estabhshment (ZFE), but the curves co- 
incide further from the source. This suggests that details about the initial 
conditions are only important only within the ZFE and that the large scale 
dynamics of the plume is determined by other parameters. 

The analysis shows a strong sensitivity to the parameter /, an increase in 
/ of 20% yielding an increase in the plume width of approximately 20%. This 
strong dependence upon I is expected and is similar to the dependence upon 
the entrainmcnt coefficient a when the entrainment hypothesis is adopted. 

The sensitivity analysis also shows a strong dependence upon the relative 
distribution width A, giving an increase in the plume width b of approximately 
40% for an increase of 20% in A. Authors using the entrainment hypothesis 
have claimed that results obtained are insensitive to the relative distribution 
width A (Bhaumik (2005), Socolofsky et. al. (2002), Wiiest (1992), Milgram 
(1983)). The discrepancy between the well established theory in the literature 
and the model presented in the present work can be related to the use of the 
entrainment hypothesis. Use of the entrainment hypothesis allows forming a 
set of non dimensional differential equations which do not depend explicitly 
on A, as shown for line sources by Ditmars and Cederwall (1974). Such a 
transformation is however not possible in the formalism based on the kinetic 
energy approach, giving an explicit dependence of A possibly explaining this 
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discrepancy. 

Even though models based on the entrainment hypothesis are relatively 
insensitive to A, little agreement is found in the literature on how this pa- 
rameter should be chosen, suggesting that the physics behind A is not clearly 
understood. The sensitivity analysis suggests that increasing the value of 
A yields an increasing plume width b and a decreasing centerline velocity 
Uc, meaning that the more the bubbles spread out, the more efficient they 
are at entraining water. This suggests the possibility that A should not be 
taken as a constant, but should depend on properties of the plume. This is 
an issue that ought to get some focus in future experimental or theoretical 
investigations. 
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Figure 5: Optimal / when A = 0.5. Solid lines represent the plume width h 
and dotted lines represent the centerline velocity Uf.- Compared with exper- 
imental data of Milgram (1983) for Q = 0.7670kg/ s. 



As an alternative to the assumption that A is dependent on flow prop- 
erties, one could investigate if the effect of choosing a different value of A 
could be incorporated in another of the free parameters of the plume. The 
sensitivity analysis suggests that a change in A could be compensated for by 
a change in I. As shown in figure [5], changing the value of A from 0.8 to 0.5, 
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changes I from 0.151 to 0.235 in order to obtain a satisfactory fit between 
theory and experiment. It is thus possible to incorporate changes in A in 
the turbulent correlation parameter / and obtain, at least qualitatively, the 
same results. If such an approach is to be used, it should be noted that the 
coefficients of equation [51] should be determined for the new values of A as 
the equation only is valid under the assumption that A = 0.8. 

6 Further issues 

1 The literature contains just a few experimental data sets that are suitable 

for comparison with model results. This is especially the case for large 
rates of release, yielding poor grounds for extrapolation of the turbu- 
lent correlation parameter I into these regions. Controlled laboratory 
experiments yield good grounds for comparison but are carried out at 
relatively small rates of release and small depths, making it difficult to 
determine the importance of dissolution. 

2 The model does not distinguish between the zone of flow establishment 

and zone of established flow. The effects of the zone of flow estab- 
lishment are believed to be handled by choosing the initial conditions 
accordingly. In the comparative study between existing models car- 
ried out by Bhaumik (2005) the approach using a densimetric Froude 
number Fr introduced by Wiiest (1992) is found to be superior to its 
alternatives. The model used in the present work shows little sensi- 
tivity to the initial Froude number Fr^ for the flow within the zone 
of established flow, suggesting that future work should focus on more 
important parameters such as the turbulent correlation parameter /, 
rather than determining exact initial conditions. 

3 In the case of strong stratification the model given in the present work will 

not be valid as effects like detrainment are not considered. This double 
plume structure can be described by means of a double plume model 
as found in McDougall (1978) and Asaeda and Imberger (1993). These 
models are based on the entrainment hypothesis and yield entrainment 
equations for each of the plumes involved, thus introducing a higher 
level of complexity. An equivalent approach for the turbulent correla- 
tion parameter / should be introduced if the kinetic energy approach 
is to advance further. 
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4 The lateral variations of the turbulent stresses are assumed to be dominat- 
ing over influence from vertical turbulence. As pointed out by Brcvik 
and Kluge (1999), the desirability of taking turbulence stresses into 
account has been emphasized by earlier workers in the field. Brevik 
and Kluge (1999) model the effects of vertical turbulence by introduc- 
ing a correction factor k. The work of these authors show that the 
factor k is interrelated to the turbulent correlation parameter /. Their 
work indicates that no definite value of k can be assigned in advance 
to an experimental situation without knowing details of the turbulence 
generating geometry of the source. A large value oi k {k — 0.3) is in- 
troduced in order to fit the model to the experimental data of Milgram 
(1983). Correlations between k and / ought to be further investigated. 

7 Concluding remarks 

In the present work, the kinetic energy approach to buoyant plumes of Bre- 
vik (1977) and Brevik and KiUie (1996) is generalized in order to allow for 
dissolution of gas. The model is compared to experiments carried out by 
Milgram (1983) and is found to reproduce experimental data with satisfac- 
tory accuracy. The results presented suggest that the model presented yields 
a good starting point for the description of the dynamics and dissolution of 
gas in the zone of established fiow for an air-bubble plume, given that the 
turbulent correlation parameter I is chosen correctly. 

The benefit of models based on the entrainment hypothesis is that the 
concept is relatively well understood and implemented for various uses. The 
major drawback of the approach is the difficulty of determining the value of 
the entrainment coefficient from parameters which are simple to measure. 

Besides being relatively simple to implement, having fast convergence 
and promising accuracy when compared to experiments, the idea of using 
a conservation equation for the kinetic energy is especially interesting from 
a physical point of view, as it provides insight into the physics driving the 
plume. The mathematical framework needed is somewhat more complicated, 
but more physics is contained in the model. The model presented should 
thus be seen as an alternative to models based on the entrainment hypoth- 
esis, with the potential of answering relevant questions without excluding 
relevant physics. Another positive aspect of the kinetic energy approach is 
that the unknown turbulent correlation parameter / can be determined from 
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parameters that are easily accessible. 

A combination of the simplicity of the entrainment hypothesis and the 
solid physical grounds of the kinetic energy approach should be sought in 
order to develop state of the art tools for future risk management of sub-sea 
gas releases, exploiting the benefits of both approaches. 
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